#include "FlowToyMCBase.h"
#include "TH1D.h"
#include <vector>
#include <iostream>
#include "TROOT.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TF1.h"
#include "TFile.h"
#include "TChain.h"
#include "TClonesArray.h"
#include "TLorentzVector.h"
#include <TStopwatch.h>

#define pi TMath::Pi()
#define PI pi
FlowToyMCBase::FlowToyMCBase() {   // Constructor

fMyRandom.SetSeed(0);

	fMass = 1;
	fYMax = 1;
	fNPar = 500;
	fVn_Input[0] = 1;
	fVn_Input[1] = 0;
	fVn_Input[2] = 0.06;
	fVn_Input[3] = 0.01;	
 	fEP_Input[0] = 0;
	fEP_Input[1] = 0;
	for (int nh=2;nh<fNH;nh++){
	fEP_Input[nh] = fMyRandom.Uniform(0,2*pi/nh); 
	}
	fdNdphi_bulk = new TF1("dNdphi_bulk", "[0]*( 1 + 2*[1]*cos(x-[4]) + 2*[2]*cos(2*(x-[5]))+2*[3]*cos(3*(x-[6])))",0,2*pi);
	fPtDist = new TF1("fPtDist", "5-x", 0, 5);  // Not realistic
 	fdNdphi_bulk->SetParameter(0, fVn_Input[0]); // underyn const : meaningless
	fdNdphi_bulk->SetParameter(1, fVn_Input[1]);  //v1  - ignore directed flow
	fdNdphi_bulk->SetParameter(2, fVn_Input[2]);  // v2 - from header file 
	fdNdphi_bulk->SetParameter(3, fVn_Input[3]);  // v3 - from header file
	fdNdphi_bulk->SetParameter(4, fEP_Input[1]); // v1 event plane  - 0
	fdNdphi_bulk->SetParameter(5, fEP_Input[2]); // v2 event plane  - random generated 
	fdNdphi_bulk->SetParameter(6, fEP_Input[3]); // v3 event palne  - random generated
	fDet_resol = 0;
	fDetName = "Default Detector";
	fArray = new TClonesArray("TLorentzVector", fNPar);

}

FlowToyMCBase::~FlowToyMCBase(){
}

void FlowToyMCBase::GeneratePhi() { 
	fPhis.clear();
	for (int ip=0; ip<fNPar; ip++) {
		double phi = fdNdphi_bulk->GetRandom();
		if (fDet_resol == 0) { fPhis.push_back(phi); }	
		else { fPhis.push_back( int(phi/fDet_resol)*fDet_resol+fDet_resol/2. ); }
	}

}
void FlowToyMCBase::Generate4Vector() { //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
		double Pt;
		double Y;
		double Phi;
		double Mt;
		double E;
		double Px;
		double Py;
		double Pz;
		FlowToyMCBase::GeneratePhi();
		for (int ipar=0;ipar<fNPar;ipar++) {
			Pt = fPtDist->GetRandom();
			Y = fMyRandom.Uniform(-1*fYMax, fYMax);
			Phi = fPhis[ipar];
			Mt = TMath::Sqrt(fMass*fMass + Pt*Pt);
			Px = Pt*TMath::Cos(Phi);
			Py = Pt*TMath::Sin(Phi);
			Pz = Mt*TMath::SinH(Y);
			E = Mt*TMath::CosH(Y);
			TLorentzVector *vt = (TLorentzVector*)fArray->ConstructedAt(ipar);
			vt->SetXYZT(Px, Py, Pz, E);			
		}
} 

void FlowToyMCBase::Init(){
	for (int nh=0;nh<fNH;nh++) {
	fQx[nh] = 0;
	fQy[nh] = 0;
	fEP_Measured[nh] = -999;
	fVn_Measured[nh] = 0;
	fVn_MeasuredReal[nh] = 0;
	fResol[nh] = 0;
	fSub2Resol[nh] = 0;
	for(int i=0;i<3;i++) {
		fSub3Resol[nh][i] = 0;
		}
	}
	if (fPhis.size() > 0) fPhis.clear();
	if (fArray->GetSize() > 0) fArray->Clear("C");
}

double FlowToyMCBase::Cal_EP(int nh) { //use this functions after Generate()
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	if (fPhis.size() == 0) {
		cout << "Generate particles first" << endl;
		exit;
	}
	fQx[nh] = 0;
	fQy[nh] = 0;
	for (int ip=0; ip<fNPar; ip++){
		fQx[nh] += TMath::Cos(nh*fPhis[ip]);
		fQy[nh] += TMath::Sin(nh*fPhis[ip]);
	}

	fEP_Measured[nh] = TMath::ATan2(fQy[nh], fQx[nh])/nh;
	return fEP_Measured[nh];
}

double FlowToyMCBase::Cal_EPSub(int nh, int ipar_from, int ipar_to) { //use this functions after Generate()
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	if (fPhis.size() == 0) {
		cout << "Generate particles first" << endl;
		exit;
	}
	if (ipar_to > fNPar) {
	cout << "Exceed Number of particle" << endl;
	exit;
	}
	fQx[nh] = 0;
	fQy[nh] = 0;
	for (int ip=ipar_from; ip<ipar_to; ip++){
		fQx[nh] += TMath::Cos(nh*fPhis[ip]);
		fQy[nh] += TMath::Sin(nh*fPhis[ip]);
	}
	double EP_MeasuredSub = TMath::ATan2(fQy[nh], fQx[nh])/nh;
	return EP_MeasuredSub;
}

double FlowToyMCBase::Cal_Vn(int nh) {
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	if (fPhis.size() == 0) {
	cout << "Generate particles first" << endl;
	exit;
	}
	if (fEP_Measured[nh] < -900) {
	cout << "EP not measured automatically call Cal_EP(nh)" << endl;
	FlowToyMCBase::Cal_EP(nh);
	}
	double Vn = 0;
	for(int ip=0; ip<fPhis.size(); ip++){
		double dphi = fEP_Measured[nh] - fPhis[ip]; 
		Vn += TMath::Cos(nh*dphi);
	}				
	fVn_Measured[nh] = Vn/fNPar;		
	return fVn_Measured[nh];
}

double FlowToyMCBase::Cal_VnReal(int nh) {
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	if (fPhis.size() == 0) {
	cout << "Generate particles first" << endl;
	exit;
	}
	double Vn = 0;
	for(int ip=0; ip<fPhis.size(); ip++){
		double dphi = fEP_Input[nh] - fPhis[ip]; 
		Vn += TMath::Cos(nh*dphi);
	}						
	fVn_MeasuredReal[nh] = Vn / fNPar;
	return fVn_MeasuredReal[nh];
}

double FlowToyMCBase::Cal_Resol(int nh){
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	if (fPhis.size() == 0) {
	cout << "Generate particles first" << endl;
	exit;
	}

	if (fEP_Measured[nh] < -900) {
	cout << "EP not measured automatically call Cal_EP(nh)" << endl;
	FlowToyMCBase::Cal_EP(nh);
	}

	double dpsi = fEP_Measured[nh] - fEP_Input[nh]; 
	fResol[nh] = TMath::Cos(nh * dpsi);
	return fResol[nh];
}
double FlowToyMCBase::Cal_Sub2Resol(int nh){
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	if (fPhis.size() == 0) {
	cout << "Generate particles first" << endl;
	exit;
	}
	double EP1 = FlowToyMCBase::Cal_EPSub(nh, 0, fNPar/2);
	double EP2 = FlowToyMCBase::Cal_EPSub(nh, fNPar/2, fNPar); 
	double dpsi = EP1 - EP2; 
	fSub2Resol[nh] = TMath::Cos(nh * dpsi);
	return fSub2Resol[nh];
}

double* FlowToyMCBase::Cal_Sub3Resol(int nh, int sub1, int sub2){
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	if (fPhis.size() == 0) {
	cout << "Generate particles first" << endl;
	exit;
	}
	double EP1 = FlowToyMCBase::Cal_EPSub(nh, 0, sub1);
	double EP2 = FlowToyMCBase::Cal_EPSub(nh, sub1, sub2); 
	double EP3 = FlowToyMCBase::Cal_EPSub(nh, sub2, fNPar);
	double dpsi1 = EP1 - EP2;
	double dpsi2 = EP2 - EP3;
	double dpsi3 = EP1 - EP3; 
	fSub3Resol[nh][0] = TMath::Cos(nh * dpsi1);
	fSub3Resol[nh][1] = TMath::Cos(nh * dpsi2);
	fSub3Resol[nh][2] = TMath::Cos(nh * dpsi3);
	return fSub3Resol[nh];
}






void FlowToyMCBase::SetNPar(int NPar){
	fNPar = (double)NPar;
}


void FlowToyMCBase::SetInputVn(double v3, double v2, double v1, double v0){
	fVn_Input[0] = v0;
	fVn_Input[1] = v1;
	fVn_Input[2] = v2;
	fVn_Input[3] = v3;
}

void FlowToyMCBase::SetEP_Measured(int nh, double EP){
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	return;
	}
	fEP_Measured[nh] = EP;	
}	



void FlowToyMCBase::SetDetResol(double Det_resol){
	fDet_resol = Det_resol;
}

void FlowToyMCBase::SetDetName(char* DetName){
	fDetName = DetName;
}

void FlowToyMCBase::SetEPSeed(int Seed){	
	fMyRandom.SetSeed(Seed);
}

void FlowToyMCBase::SetEP() {
 	fEP_Input[0] = 0;
	fEP_Input[1] = 0;
	fEP_Input[2] = fMyRandom.Uniform(0,2*pi/2); 
	fEP_Input[3] = fMyRandom.Uniform(0,2*pi/3);
	fdNdphi_bulk->SetParameter(4, fEP_Input[1]); // v1 event plane  - 0
	fdNdphi_bulk->SetParameter(5, fEP_Input[2]); // v2 event plane  - random generated 
	fdNdphi_bulk->SetParameter(6, fEP_Input[3]); // v3 event palne  - random generated
}

void FlowToyMCBase::SetPtDist(TF1* PtDist) {   
	fPtDist = PtDist;
}

void FlowToyMCBase::SetRapidityRange(double YMax) {
	fYMax = YMax;
}
  
void FlowToyMCBase::SetMass(double Mass) {
	fMass = Mass;
}

double FlowToyMCBase::GetQx(int nh){
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	return 0;
	}
	return fQx[nh];
}

double FlowToyMCBase::GetQy(int nh){
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	return 0;
	}
	return fQy[nh];
}


double FlowToyMCBase::GetEP_Input(int nh){
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	return fEP_Input[nh];
}

double FlowToyMCBase::GetEP_Measured(int nh){
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	return fEP_Measured[nh];
}

double FlowToyMCBase::GetVn_Input(int nh){
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	return fVn_Input[nh];
}

double FlowToyMCBase::GetVn_Measured(int nh){
	if (nh >= fNH) {
	cout << "Exceed nh range" << endl;
	exit;
	}
	return fVn_Measured[nh];
}

double FlowToyMCBase::GetPhi(int iPar){
	return fPhis[iPar];
}

double FlowToyMCBase::GetPx(int iPar) {
	TLorentzVector *vt = (TLorentzVector*)fArray->At(iPar);
	double Px = vt->X();
	return Px;
}

double FlowToyMCBase::GetPy(int iPar) {
	TLorentzVector *vt = (TLorentzVector*)fArray->At(iPar);
	double Py = vt->Y();
	return Py;
}

double FlowToyMCBase::GetPz(int iPar) {
	TLorentzVector *vt = (TLorentzVector*)fArray->At(iPar);
	double Pz = vt->Z();
	return Pz;
}

double FlowToyMCBase::GetE(int iPar) {
	TLorentzVector *vt = (TLorentzVector*)fArray->At(iPar);
	double E = vt->T();
	return E;
}

double FlowToyMCBase::GetPt(int iPar) {
	TLorentzVector *vt = (TLorentzVector*)fArray->At(iPar);
	double Pt = vt->Pt();
	return Pt;
}

double FlowToyMCBase::GetRapidity(int iPar) {
	TLorentzVector *vt = (TLorentzVector*)fArray->At(iPar);
	double Y = vt->Rapidity();
	return Y;
}

TClonesArray* FlowToyMCBase::GetParticleArray() {
	return fArray;
}


void FlowToyMCBase::Print(){
 cout << "what?" << endl;
}

double FlowToyMCBase::GetNPar(){
	return fNPar;
}


